/*
 * 2-dimensional Ising model, with Markov chain in Monte Carlo method
 * 	H = - J \sum_{ij} sz_i sz_j - muB \sum_i sz_i,
 * where J is the interaction strength between two atoms, negative/positive for ferromagnetic/anti-ferromagnetic material, k is the Boltzman constant, T is the temperature, we set J as unit.
 */

#include<iostream>
#include<stdlib.h>
#include<time.h> 
#include<fstream>
#include<cmath>
using namespace std;

void Ising3d(int length, int L, double kT, double muB, double & sz_ave, double & e_ave, double & e2_ave){

	long i,j,k, Sz, S_neighbour, count=0;
	double J=1, E, dE, sz_sum=0, e_sum=0, e2_sum=0;
	int *spin=new int [L*L*L];//spin of each atom, which can only be 1 or -1
	bool flag_flip;

	for(i=0;i<L*L*L;i++) //initially the lattice is magnetizied.
		spin[i]=1;

	srand((unsigned)time(0));//initialize rand() with a random seed, according to the real time when this code runs.
	Sz= L*L*L;//initial value of total spin
	E=-J*3*L*L*L - muB * Sz;//initial value of total energy of the system

	for(j=0;j<length;j++){

		for(i=0;i<L*L*L;i++){

			S_neighbour = spin[ i%L ==0 ? i+L-1 : i-1 ] 
				+ spin[ (i+1)%L==0 ? i-(L-1) : i+1 ]
			       	+ spin[ i%(L*L) >= (L-1)*L ? i-(L-1)*L : i+L ] 
				+ spin[ i%(L*L) < L ? i+(L-1)*L : i-L ]
				+ spin[ i<(L*L) ? i+(L-1)*L*L : i-L*L ]
				+ spin[ i>=(L-1)*L*L ? i-(L-1)*L*L : i+L*L ];
			dE= 2 * J * S_neighbour * spin[i] + 2 * muB * spin[i];

			flag_flip=false;
			if(dE<0){
			       flag_flip=true;
			}
			else{
				if( (double)rand()/RAND_MAX < exp(-dE/kT) ){// rand()/RAND_MAX uniformly distributes in [0,1], so it has possibility exp(-dE/kT) to be less than exp(-dE/kT)
					flag_flip=true;
				}
			}
		
			if(flag_flip){
				E += dE;
				spin[i] *= -1;//flip the spin
				Sz += 2*spin[i];// total spin and total energy changes
			}
		}
		// after one round, collects data

		if(j>0.1*length){// We assume the system is in equilibrium after 0.1*length
			sz_sum += (double)Sz/(L*L*L);
			e_sum += E/(L*L*L);
			e2_sum += E*E/(L*L*L);

			count ++;
		}
	}

	sz_ave = sz_sum/count;
	e_ave = e_sum/count;
	e2_ave = e2_sum/count;

	delete [] spin;
}

int main(){

	int length = 1E3, L=8;
	double sz_ave, e_ave, e2_ave, muB=0;

	ofstream fp("Ising3d.txt");
	for(double kT=0.1; kT<=7; kT+=0.1){
		Ising3d(length, L, kT, muB, sz_ave, e_ave, e2_ave);
		fp<< kT <<"\t"<< sz_ave << "\t" << e_ave <<endl;
	}
	fp.close();

	return 0;
}
